Acetylcholinesterase (AChE) contains a nucleophilic serine that reacts with acetate esters and organophosphorus compounds - phosphates and phosphonates. Organophosphonates (OP) exist in two enantiomeric forms. The high activity of them with AChE is most probably caused by the fact that they are transition-state tetrahedral analogs of substrate hydrolysis, so ground-state tetrahedral symmetry of these agents would be expected to determine reactivity. On the other hand, if the 3-d geometries of these agents compliment the active center of AchE surface, then the nature of groups surrounding phosphorus would determine reactivity. To find out the answer to that question one can model the process of inhibition of OP using enantiomeric pairs of cycloheptyl-, isopropyl and 3,3-dimethylbutyl methylphosphothionates. A number of experimental works were done to study interactions of these OP with AChE. Several calculations were done to study the mechanism of OP - AchE (Buche) interacti ons. A successful study of OP binding process to AchE requires knowledge of the reliable initial conformations of all OP involved in this process. So far that conformations were constructed using MOPAC optimizer from starting conformation created by InsighII Builder program. But MOPAC approach does not give the opportunity to calculate the structure of OP in solution and does not give the most precise results of optimization. To avoid that problem one can use the DMol/COSMO program which combines accurate electronic structure calculations using DFT method with simulation of the solvent effects via COSMO model. In this model, the solute molecule is embedded into a cavity surrounded by the solvent, represented as a dielectric continuum. The polarization of the dielectric continuum by the charge distribution of the solute results in a charge distribution on the cavity surface. In the COSMO method, the surface charges are obtained directly from the electrostatic potential on the cavity surface. Since the DMol/COSMO energy is fully variational, accurate gradients with respect to the solute coordinates can be calculated and therefore one can perform optimization of the molecular structure in solvent environment. As a first step, the optimized structures of 6 organophosphonates pairs were calculated using DMol/COSMO program in gas phase and in water. Relative energies of conformers were found and the charge distribution due to the solvent effects was calculated. Further calculations are planned to estimate the transition states of such OP reactions with active serine of AchE and other esterases.